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Abstract 

Holographic  Particle  Image  Velocimetry  (HPIV)  is  a  novel  technique  for  measuring 
the  complete  fluid  flow  around  a  body.  Advances  in  computing  power  make  this 
technique  practicable  for  the  first  time. 

Currently  popular  techniques  for  experimentally  determining  fluid  flow  around 
a  test  body  rely  on  measuring  the  flow  at  a  single  point  and  moving  the  sample 
point  during  the  experiment.  This  implicitly  time  averages  the  data.  HPIV  can 
measure  the  flow  field  in  a  volume  nearly  instantaneously  by  sequentially  "holo- 
graphing" the  seeded  flow  field. 

Until  recently  the  only  practical  method  of  processing  the  HPIV  data  was  to 
physically  reconstruct  the  optical  image  from  the  hologram  and  scan  through  the 
resulting  volume  using  microscope  optics.  This  is  a  very  slow  process. 

Recent  advances  in  computing  power  make  it  possible  to  directly  analyze  the 
digitized  holograph  (digitized  using  a  high  resolution  scanner)  This  reduces  the 
post-processing  time  from  several  hours  to  less  than  1-hour  with  much  less  expen- 
sive equipment. 

The  primary  contribution  of  this  thesis  is  an  advanced  multi-platform  graphical 
user  interface  hologram  synthesizer  for  use  in  calibrating  digital  HPIV  reconstruc- 
tion algorithms. 

Thesis  Supervisor:  Jerome  H.  Milgram 
Title:  Professor  of  Ocean  Engineering 
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Chapter  1 


Introduction 


1.1  Overview 

This  thesis  presents  an  advanced  multiplatform  graphical  user  interface  hologram 
synthesizer  to  produce  calibration  holograms  of  seeded  flow  fields.  The  synthe- 
sizer, HoloSynth,  runs  under  X11R6  with  OpenGL,  and  Microsoft  Windows  NT. 
Complete  source  code  is  contained  in  Appendix  A.  A  brief  introduction  to  parti- 
cle holography  and  holographic  particle  velocimetry  is  presented  as  well  as  results 
from  computational  reconstruction  of  the  hologram  generated  by  the  synthesizer. 

1.2  Motivation 

Until  very  recently,  particle  image  velocimetry  (PIV)  was  limited  to  measuring 
fluid  flow  velocities  only  in  a  single  plane.  A  sheet  of  laser  light  was  passed 
through  the  flow  which  contained  neutrally  buoyant  particles.  A  double  exposure 
of  the  particles  onto  either  film  or  a  Charge  Coupled  Device  (CCD)  camera  allows 
the  two  velocity  components  in  the  illuminated  plane  to  be  determined. 

The  PIV  method  has  recently  been  extended  to  three  dimensions  by  Katz  et.  al. 
[1]  by  making  a  hologram  of  the  flow  field.  Three  dimensional  Holographic  Particle 
Image  Velocimetry  (HPIV)  uses  either  in-line  or  off-axis  holographic  techniques 
(see  Section  2.2).  Once  the  hologram  has  been  recorded,  the  image  is  reconstructed 
by  reversing  the  recording  process  and  using  a  three  dimensional  optical  scanner  to 


1.3.     REQUIREMENTS  10 

digitize  the  image  for  analysis  by  computer.  This  process  can  take  several  hours. 

An  alternative  to  scanning  the  reconstructed  image  is  to  digitize  the  hologram 
and  allow  a  computer  to  reconstruct  the  image  in  software.  The  recent  advances  in 
computing  power  should  make  this  process  far  faster  than  the  mechanical/optical 
scanning  process  now  used. 

Digitizing  the  hologram  can  be  done  in  two  ways.  The  initial  method  avoids 
the  three  dimensional  scanning  of  the  reconstructed  image  by  using  a  high  perfor- 
mance flat  bed  scanner  to  digitize  an  enlargement  of  the  transparency  produced 
in  the  normal  holographic  process.  Alternately,  a  high  resolution  scanner  can  be 
used  directly  on  the  hologram  transparency.  The  flat  bed  scanning  process  is  well 
understood,  since  the  transparency  can  be  enlarged  to  overcome  resolution  prob- 
lems with  the  scanner.  However,  it  still  requires  recording  the  hologram  on  film 
with  the  attendant  film  processing. 

A  more  advanced  technique  would  be  to  capture  the  hologram  using  a  Charge 
Coupled  Device  (CCD)  Camera.  The  hologram  would  be  immediately  available 
in  digital  form,  and  no  film  or  processing  would  be  involved. 

Regardless  of  the  hologram  digitization  technique  the  computational  machinery 
for  reconstructing  the  image  must  be  calibrated.  It  far  easier  to  produce  a  synthetic 
hologram  for  calibration  than  to  set  up  a  tightly  controlled  experiment  and  actually 
record  a  hologram. 

1.3     Requirements 

The  initial  goals  of  computationally  reconstructed  HPIV  are  to  be  able  to  accu- 
rately determine  the  flow  field  in  a  cubic  volume  a  few  inches  on  a  side.  The  fluid 
speeds  will  be  typical  of  those  used  in  hydrodynamics  research.  The  quantifiable 
requirements  are: 

•  Must  be  able  to  resolve  a  10/Lmi  particle.  This  corresponds  to  the  most 
conveniently  sized  neutrally  buoyant  seed  particles  used  for  Laser  Doppler 
Anemometry. 
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•  Accurate  coverage  within  a  125cm3  (5cm  cube)  volume.  This  will  rise  as  the 
technique  is  perfected,  and  computing  power  increases. 

•  Unambiguously  detect  fluid  speeds  <  lOm/s.  This  accounts  for  free  stream 
plus  allowance  for  turbulence  and  model  scale  factors. 

The  first  and  second  requirements  dictate  the  size  and  resolution  of  the  recording 
device. 

Since  holography  is  essentially  a  linear  process,  multiple  exposures  can  be  used 
instead  of  a  rapid  succession  of  exposures.  So  rather  than  design  a  system  to 
rapidly  capture  two  separate  holograms,  the  recording  medium  is  exposed  two 
or  more  times  without  being  altered  between  exposures.  This  results  in  multiple 
real  images  for  each  particle  in  the  volume,  the  analysis  process  must  identify  the 
multiple  images  with  a  particular  particle.  This  alleviates  the  need  for  a  high 
speed  film  traversal  mechanism  and  removes  CCD  output  data  rate  as  a  concern. 

The  third  requirement  sets  the  pulse  lengths  and  pulse  intervals  of  the  laser 
used  to  generate  the  hologram.  If  the  pulses  are  too  long  the  diffraction  patterns 
will  be  smeared  from  the  particle  motion.  The  greater  the  interval  between  pulses, 
the  harder  it  will  be  to  uniquely  identify  particle  images  with  a  single  particle.  The 
pulse  length  and  interval  requirements  determine  the  recording  device  sensitivity 
and  laser  output  power. 
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Chapter  2 


Review  of  Holographic  Technique 


2.1      Overview 

A  hologram  is  a  two  dimensional  recording  that  can  produce  a  three  dimensional 
image.  Typically  electromagnetic  (light)  waves  are  recorded,  but  recording  of 
acoustic  waves  and  matter  waves  have  also  been  made[2].  Electromagnetic  holog- 
raphy differs  from  normal  photography  in  that  it  encodes  phase  information  as 
well  as  intensity  information  on  the  recording  medium. 

Recording  media  respond  only  to  the  intensity  of  the  incoming  wavefront.  not 
its  phase.  To  encode  the  phase  using  normal  media,  holography  depends  on  using  a 
coherent  reference  beam  to  interfere  with  a  scattered  or  diffracted  wave,  or  "object 
wave"',  from  the  objects  under  observation.  As  a  result  of  this  interference,  the 
intensity  at  the  recording  plane  depends  on  both  the  phase  and  amplitude  of  the 
object  wave. 

The  original  image  is  reproduced  simply  by  illuminating  the  hologram  with 
coherent  light  using  geometry  similar  to  that  used  to  record  the  hologram.  If  the 
wavelength  used  for  reconstruction  differs  from  the  wavelength  used  for  recording, 
the  image  will  be  enlarged  or  reduced. 

The  reconstruction  process,  while  simple,  is  not  amenable  to  rapid  quantitative 
data  collection.  Quantitative  analysis  generally  involves  using  three  dimensional 
optical  scanning  of  the  reconstructed  image.    The  ability  to  computationally  re- 
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Figure  2-1:  Basic  in-line  hologram  geometry 


construct  the  image  would  dramatically  increase  the  utility  of  holographic  mea- 
surement techniques. 


2.2     In-line  (Gabor)  Holography 

In-line  holography,  or  Gabor  Holography  after  its  inventor,  describes  a  holographic- 
apparatus  in  which  the  reference  beam  is  usually  oriented  perpendicular  to  the 
recording  media  and  illuminates  the  objects  from  behind.  It  is  not  a  useful  tech- 
nique for  objects  that  are  large  with  respect  to  the  recording  medium,  but  if  one 
is  interested  primarily  in  the  three  dimensional  position  of  small  objects  and  not 
their  three  dimensional  shape,  in-line  holography  can  be  very  effective.  It  is  com- 
monly used  to  size  particulates  and  aerosols  suspended  in  fluids,  and  it  has  been 
used  for  traditionally  reconstructed  HPIV[3]. 

Figure  2-1  shows  the  basic  geometry  and  coordinate  system  for  recording  a 
plane  in-line  hologram.  The  object  is  illuminated  by  the  reference  source.  The 
diffracted  wave,  or  object  wave,  interferes  with  the  wave  from  the  reference  source, 
or  reference  wave.  This  interference  pattern  is  recorded  by  the  medium  in  the 
recording  plane  and  constitutes  the  hologram  itself. 

If  the  reference  and  object  complex  wave  amplitudes  measured  at  the  recording 
plane  are  given  by  R  (x,  y)  and  O  (x,  y)  respectively,  the  intensity  at  the  recording 
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plane  is  given  by. 


I(x,y)  =  \R(x,y)  +  0(x,y)f 


(2.i: 


Since  the  recording  plane  is  perpendicular  to  the  plane  reference  wave,  the 
complex  amplitude  of  the  reference  wave  is  a  constant,  r.  Expanding  (2.1)  yields 


I{x,y)    =     \r  +  0{x,y)\ 


=    r2  +  \0{x,y)\"  +rO(x,y)  +  rO*(x,y) 


(2.2; 


The  intensity  is  recorded  in  any  acceptable  manner,  usually  as  a  positive  trans- 
parency on  film.  Most  reasonable  recording  media  respond  nearly  linearly  to  in- 
cident wave  intensity.  The  transmittance,  t(x,y),  of  a  positive  transparency  can 
be  adequately  described  by 


t{x,y)    =    rb-PTI{x,y) 


=    n-pr 


r2  +  \0(x,y)\~  +  rd(x,y)+r6*(x,y)  (2.3) 


where  rb  is  a  constant  background  transmittance,  T  is  the  exposure  time,  and  (3 
is  a  constant  of  the  recording  medium. 

Figure  2-2  shows  an  in-line  hologram  being  physically  reconstructed.     The 
wavefront  transmitted  by  the  hologram  is  simply  the  incident  wave  multiplied  by 
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the  transmittance  of  the  hologram  and  can  be  written, 

ijj{x,y)    =    rt{x,y) 

=    r(n  +  PTr2)+(3Tr\d{x,y)\2 

+pTr20  (x,  y)  +  0Tr2O*  (x,  y)  (2.4) 

The  first  term  in  (2.4)  represents  the  uniformly  attenuated  reference  wave.  The 
second  term  is  normally  very  small,  since  typically  O  (x,  y)  <C  r.  It  is  called  the 
"halo''  and  under  proper  conditions  can  be  neglected.  The  third  term  is  identical 
to,  within  a  multiplicative  constant,  with  the  original  scattered  wave.  It  forms  an 
image  z0  from  the  hologram  on  the  side  opposite  the  observer.  This  is  the  virtual 
image,  which  is  what  the  observer  perceives.  The  fourth  term  represents  a  wave- 
front  similar  to  the  originally  scattered  wavefront,  but  with  opposite  curvature.  It 
converges  to  form  a  real  image  z0  from  the  hologram,  between  the  hologram  and 
the  observer.  [4] 

The  observer  then  sees  the  virtual  image,  with  a  superimposed,  out-of-focus 
real  image  along  with  the  diminished  reference  beam.  The  two  images  instead  of 
one  and  are  the  biggest  drawbacks  to  in-line  holography. 

2.3      Off-axis  (Leith-Upatneiks)  Holography 

A  simple  way  to  separate  the  out-of-focus  real  image  from  the  virtual  image  is  to 
provide  a  separate  path  for  the  reference  beam  and  the  object  wave.  Figure  2-3 
shows  the  basic  geometry  for  recording  an  off- axis  hologram. 

Because  the  reference  beam  is  no  longer  perpendicular  to  the  recording  plane, 
it  cannot  be  simplified  to  a  constant,  rather,  it  can  be  expressed  as, 

R  (x,  y)  =  r  exp(ji'27r^rx)  (2.5) 

where  £r  =  (sin#)/A,  since  only  the  phase  of  the  reference  beam  varies  across  the 
recording  plane. 
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Figure  2-3:  Off-axis  hologram  geometry 

In  a  fashion  similar  to  the  derivation  of  (2.2)  the  intensity  at  the  recording 
plane  is  found  to  be. 

/  (:r,  y)  =  r2  \0  (x,y)\2  +  2r\0(x,y)\  cos  [2^rar  +  4>{x,y)]  (2.6) 

Figure  2-4  shows  the  geometry  for  an  off-axis  hologram  reconstruction.  As- 
suming the  recording  medium  is  linear,  the  transmitted  wave  complex  amplitude 
is  found  in  a  similar  manner, 

■0  (x,  y)    =    rbr  exp(j2n£rx) 

+PTr  \0  {x,y)\   exp(j2n£rx) 

+!3Tr20(x,y) 

+[3Tr2d{x,y)exp(j47r£rx)  (2.7) 

Once  again,  the  first  term  in  (2.7)  is  simply  the  attenuated  reference  beam. 
The  second  term  is  a  "halo"  term,  similar  to  the  neglected  term  in  (2.7).  The 
third  term  is  identical  to  the  original  scattered  wave  and  forms  a  virtual  image 
in  the  same  position  as  the  original  object.  The  fourth  term  is,  again,  the  real 
image,  this  time  shifted  from  the  optical  axis  by  approximately  twice  the  angle 
the  reference  beam  makes  with  the  recording  plane.  Hence  the  out  of  focus  real 
image  does  not  interfere  with  the  vitual  image. 
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Figure  2-4:  Off-axis  hologram  reconstruction  geometry 

At  first  blush  it  would  seem  off-axis  holography  is  the  answer  to  the  virtual 
image  problem  for  HPIV.  If  purely  physical  reconstruction  is  to  be  used,  then  off- 
axis  holography  is  eminently  practicable.  The  fatal  flaw  comes  when  considering 
digitizing  the  hologram  for  use  in  computational  reconstruction.  The  reference 
beam  now  introduces  a  rapidly  fluctuating  intensity  term  that  quadruples  the 
resolution  requirement  of  the  recording  medium  and  digitization  system.  For  this 
reason  off-axis  holography  will  not  be  discussed  further  in  this  thesis. 

2.4     Diffraction 

Holography  is  fundamentally  about  diffraction.    Goodman [5]  provides  a  detailed 
derivation  of  the  Huygens-Fresnel  principle  for  a  wave  from  constant  z, 

z     ff,i.r„  „AexP(j^roi) 

„2 

01 


t/j(xh,yh)  =  —  J  J  tf>x(x,y) —2 dxdy  (2.8) 


>: 


where,  r0i  =  yj  z2  +  (xh  —  x)2  +  (yh  —  y)2,  and  Xf,  and  ;<//,  are  coordinates  in  the 
recording  plane. 

This  expression,  while  physically  meaningful,  is  not  easily  evaluated.  Two 
fundamental  simplifications  are  common  depending  on  the  geometry  of  the  optical 
problem  being  analyzed.  The  approximations  are  the  Fresnel  and  the  Fraunhofer 
approximations,  valid  in  the  near  and  far  fields,  respectively. 
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2.4.1      The  Fresnel  Approximation 

A  simpler  and  more  useful  expression  of  the  Huygens- Fresnel  integral  can  be  found 
by  simplifying  r0i  using  the  paraxial  approximation,  z2  ^>  x2  +  y2. 
If /;  <  1.  then, 

xfTTb=l  +  h-lb2  +  ---  (2.9) 

2  o 


Using  (2.9)  with  b  =  f^)2  +  (**f*) 


roi     =     V  -2  +  (xh  ~  x)2  +  {yh  -  y)2 


(2.10) 


Substitution  of  (2.10)  into  (2.8)  yields  little  simplification  unless  carefully  cho- 
sen terms  are  dropped.  For  the  r^  term  in  the  denominator  of  (2.8)  it  is  considered 
acceptable  to  drop  all  terms  but  z.  However,  small  changes  in  the  r0i  occurring 
in  the  numerator  can  cause  significant  changes  in  the  value  of  the  exponential 
and  cannot  be  omitted.  These  large  changes  are  generated  in  part  by  the  mul- 
tiplication of  r0i  by  A:,  usually  a  very  large  quantity,  and  the  sensitivity  of  the 
exponential  to  phase  changes  in  its  argument.  The  resulting  simplification  yields 
the  Fresnel- Kirchoff  Diffraction  Integral: 


oo 


■&(xh,yh)     =     (^U{x,y)explj—[{xh-xj2  +  {yh-y)2]ldxdi2.n) 

—  oo 

oc 

=     tJ-^M+vl)    f  f  U{x,y)eJ^(x2+y2)e-J%{XhX+yhy)dxdy(2.l2) 


Where  (2.12)  differs  from  (2.11)  only  through  factoring  the  term  exp  [jj  (x2h  +  yfj] 
outside  the  integral.   This  approximation  holds  when  the  observer  is  in  the  near 
field,  or  Fresnel  zone,  of  the  diffracting  aperture.    Typical  HPIV  geometry  will 
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place  the  recording  medium  in  the  far  field  of  the  individual  particles  but  in  the 
near  field  of  the  ensemble. 

2.4.2      The  Fraunhofer  Approximation 

Additional  simplifications  to  the  Fresnel  Diffraction  Integral  can  be  made  if  the 
system  meets  the  Fraunhofer  approximation. 


fc(*2+S/2) 


max 


z>  — -_^l  (2.13; 

Simplifying  (2.12)  yields, 


iP{xh,yh)  =  e—e3U*l+yl)  /I  U{x,y)e-J^(Xh:r+yhy}dxdy  (2.14) 


which  is  valid  in  the  far  field.   The  integral  is  recognizable  as  a  two  dimensional 
Fourier  transform. 

2.4.3      Convolution  interpretation  of  Fresnel-Kirchoff  Inte- 
gral 

The  Fresnel-Kirchoff  Integral, 

oo 

iPixhiVh)  =  e—,         U(x,y)explj—  [{xh  -  xf  +  (yh  -  yf]  \  dx  dy       (2.15) 
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can  easily  be  viewed  as  a  convolution  with  kernel, 

(2.1G) 
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yielding 
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ipz{xh,yh)  =    //  U{x,y)hz(xh  -x,yh  -  y)dx dy  =  U{xh,yh)  0  hz{x,y)      (2.17) 
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This  expression  determines  the  field  diffracted  from  a  planar  object,  or  an  ensemble1 
of  co-planar  planar  objects.  If  a  large  number  of  very  small  non-coplanar  objects 
are  in  the  field  of  view,  (2.17)  can  be  extended: 

ip(xh,yh)  =  ^2  Un(xh,yh)  ®  hZn(xh,yh)  (2.18) 

where  the  sum  is  taken  over  n  planes  within  the  volume.  Each  object  plane. 
Un{x,  y),  contains  at  least  one  diffracting  object  and  is  located  zn  from  the  record- 
ing plane.  The  planes  need  not  be  evenly  spaced. 

2.5      Laser  Characteristics 

2.5.1      Fundamentals 

At  the  atomic  level,  light  is  always  generated  by  an  atom  transitioning  from  one 
energy  state,  E0,  to  a  lower  energy  state,  E\ .  When  this  transition  occurs  the  atom 
emits  a  photon  with  frequency  v  =  (EQ  —  E\)/h,  where  h  is  Planck's  constant. 
This  emission  is  a  probabilistic  occurrence  and  may  occur  spontaneously  provided 
the  atom  has  sufficient  internal  energy.  This  is  spontaneous  emission. 

A  second  form  of  emission,  first  suggested  by  Einstein  in  1917,  is  called  stim- 
ulated emission.  In  the  case  of  stimulated  emission,  a  photon  collides  with  an 
excited  atom.  The  incident  photon  must  have  exactly  the  same  energy  as  the 
transition  the  atom  will  undergo.  When  this  condition  is  met,  the  atom  will  make 
the  transition,  emitting  a  second  photon.  The  emitted  photon  is  in  phase  with 
the  incident  photon  and  travels  in  the  same  direction. 

Stimulated  emission  is  fundamental  to  laser  operation.  A  material  is  induced  to 
have  a  majority  of  its  atoms  in  an  excited  state  (called  a  "population  inversion"). 
Since  the  majority  of  the  atoms  in  the  material  are  excited  to  the  same  state, 
they  have  identical  transition  energies.  Initial  spontaneous  emission  by  a  small 
number  of  the  atoms  will  stimulate  emission  in  others.  The  stimulated  photons 
will  in  turn  stimulate  further  emission  in  a  run-away  chain  reaction.  This  in  itself 
will  not  cause  light  amplification,  but  if  the  material  is  mounted  within  an  optical 
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cavity  (two  parallel  reflectors),  the  stimulated  emission  will  preferentially  orient 
itself  along  the  axis  of  the  cavity.  Radiation  not  parallel  to  the  axis  will  be  lost 
to  the  surroundings  and  not  contribute  to  the  electromagnetic  wave  building  up 
in  the  cavity.  Since  the  stimulated  emission  is  always  in  phase,  the  radiation 
traveling  along  the  cavity  axis  constructively  reinforces  and  the  amplitude  rapidly 
increases  forming  a  longitudinal  vibration  mode  within  the  cavity.  If  an  end  of 
the  cavity  is  only  partially  reflective,  then  some  of  the  coherent  amplified  signal 
will  be  released,  this  is  the  usable  laser  beam. 

2.5.2      Establishing  the  Population  Inversion 

The  process  of  establishing  the  population  inversion  is  called  'pumping'.  Pumping 
can  be  accomplished  using  flash  lamps,  usual  in  the  case  of  ruby  lasers,  or  electrical 
discharge,  as  in  the  case  of  most  gas  lasers.  Photochemical  reactions  can  also 
stimulate  the  population  inversion  but  that  technique  is  very  rare.  [6] 

The  cavity  material  characteristics  fundamentally  determine  what  type  of  pump- 
ing will  be  appropriate.  Continuous  beam  lasers  must  obviously  use  a  sustainable 
pumping  technique,  such  as  tungsten  filament  incandescent  lamps,  or  the  afore- 
mentioned electrical  discharge. 

The  length  of  the  pulse  in  non-switched  pulse  lasers  is  directly  related  to  the 
length  of  the  pulse  length  of  the  pump.  Flash  lamp  pumped  ruby  lasers  typically 
provide  pulses  on  the  order  of  several  milliseconds  long.  The  length  of  the  pulse 
also  varies  from  pulse  to  pulse  as  the  condition  of  the  flash  lamp  changes. 

The  minimum  time  between  lasers  pulses  in  non-switched  lasers  is  directly 
related  to  the  requirements  of  the  pump.  Flash  lamps  typically  must  cool  imme- 
diately after  use.  The  amount  of  cooling  is  directly  related  to  the  power  output 
of  the  lamp,  and  can  be  on  the  order  of  minutes.  This  limits  the  pulse  repetition 
frequency  of  the  unswitched  laser. 

Switching  can  reduce  pulse  lengths  to  femtoseconds  (lxlO-15  seconds)  and 
pulse  intervals  similarly.  Switching  is  discussed  in  some  detail  in  Section  2.5.5. 
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2.5.3  Energy  Levels 

When  common  laser  materials  are  pumped,  the  resulting  atomic  population  usu- 
ally has  many  more  than  one  possible  transition  state  available.  The  likelihood 
that  an  atom  undergoes  a  particular  transition  depends  on  the  quantum  mechan- 
ics of  the  particular  material.  In  any  event,  when  the  population  inversion  starts 
to  break  down  there  can  be  a  large  number  of  longitudinal  oscillation  modes, 
each  producing  a  separate  wavelength  in  the  cavity.  Typically,  a  small  number  of 
closely  spaced  wavelengths  will  predominate,  but  the  resultant  beam  is  not  strictly 
monochromatic  without  special  modifications  to  the  laser  (see  Section  2.5.7).  [6] 

2.5.4  Coherence  Length 

The  preceding  discussion  leaves  the  impression  that  all  light  generated  in  the  laser 
cavity  is  coherent.  This  is  not  the  case.  Not  all  stimulated  photons  are  descended 
from  the  same  spontaneously  emitted  photon,  and  the  Heisenberg  uncertainty 
principle  places  a  lower  limit  on  the  precision  with  which  the  wavelengths  are 
reproduced.  So  even  a  laser  oscillating  in  only  a  single  longitudinal  mode  has 
some  degree  of  incoherence  in  its  output. 

One  measure  of  the  quality  of  a  laser's  output  is  termed  "coherence  length'' . 
This  length  is  the  maximum  beam  path  length  within  which  the  beam  can  coher- 
ently interfere  with  itself.  The  practical  ramification  of  this  is  that  the  holographic- 
system  must  have  all  path  lengths  from  laser  output  to  recording  planes  less  than 
the  coherence  length  of  the  laser.   [7] 

2.5.5  Switching 

It  is  clear  that  the  amplification  properties  of  the  laser  cavity  are  critically  depen- 
dent on  the  alignment  and  quality  of  the  end  mirrors.  Misalignment  or  surface 
imperfections  can  destroy  the  amplification  effect  and  turn  the  laser  into  an  ex- 
pensive flash  bulb. 

This  sensitivity  can  be  turned  to  the  laser  users  advantage.  If  the  mirror  prop- 
erties can  be  accurately  controlled  then  the  amplification  process  can  be  controlled 
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independently  of  the  pumping  process. 

Borrowing  a  term  from  electrical  engineering,  the  ability  of  the  cavity  to  res- 
onate is  called  it's  'Q'.  Systems  that  switch  the  amplification  process  on  and  off 
by  controlling  the  resonance  qualities  of  the  cavity  are  called  'Q-switched'  or  'Q- 
spoiled"  lasers. 

Q-switching  has  a  subtle  benefit.  The  population  inversion  does  not  occur 
instantaneously  at  the  beginning  of  the  pumping  sequence.  If  the  cavity  cannot 
oscillate  due  to  low  Q,  the  population  inversion  can  build  to  a  much  higher  level 
than  required  for  lasing.  If  the  Q  is  suddenly  increased,  allowing  the  cavity  to 
oscillate,  the  peak  energy  in  the  pulse  can  be  much  higher  than  if  the  cavity  was 
allowed  to  oscillate  continuously  through  the  pumping  sequence.  So  the  power 
available  in  a  Q-switched  pulse  can  be  much  higher  than  in  an  unswitched  pulse.  [6] 

An  additional  benefit  of  Q-switching  is  the  ability  to  produce  multiple  pulses 
out  of  a  single  pumping  process.  This  allows  accurate  control  of  continuously 
pumped  lasers,  and  multiple  rapid  pulses  out  of  short  pulse  pumped  lasers  such 
as  ruby  lasers. 

2.5.5.1      Electrooptical  Q-Switching 

Several  physical  phenomena  allow  the  electrical  modulation  of  an  optical  signal. 
The  most  commonly  used  in  Q-switching  lasers  is  embodied  in  the  Pockels  Cell. 

Pockels  Cells  are  formed  from  a  crystal,  typically  potassium  dihydrogen  phos- 
phate (KDP).  KDP  has  two  distinct  optical  axes.  With  an  electrical  voltage 
applied  across  the  crystal  the  index  of  refraction  on  one  axis  is  different  than  the 
other.  An  incident  polarized  wave  is  split  into  two  portions,  one  portion  traveling 
more  slowly  than  the  other.  The  two  waves  are  recombined  upon  exiting  the  cell, 
resulting  in  a  wave  that  is  elliptically,  circularly,  or  linearly  polarized  depending 
on  the  applied  voltage. 

The  two  voltages  of  interest  for  Pockels  Cell  Q-switching  are  termed  the  quarter- 
wavelength  (A/4)  voltage,  and  the  half- wavelength  (A/2)  voltage.  The  quarter- 
wavelength  voltage  will  result  in  a  circularly  polarized  beam  at  the  output  of  the 
cell.    The  half- wavelength  voltage  will  result  in  a  wave  linearly  polarized  90°  to 
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Figure  2-5:  Pockels  Cell  Q-switch  at  (a)  quarter-wave  and  (b)  half  wave  retardation 
voltage. 


the  input  wave. 

The  quarter-wavelength  retardation  system  shown  in  Figure  2-5(a)  is  a  'pulse 
off'  Q-switch.  During  the  initial  phase  of  pumping,  the  A/4  voltage  is  applied  to 
the  cell.  The  output  of  the  Pockels  cell  is  circularly  polarized,  after  reflection  the 
circularly  polarized  beam  reenters  the  Pockels  cell  and  emerges  linearly  polarized 
90°  to  the  input  beam.  The  polarizer  then  prevents  the  beam  from  reentering 
the  laser  rod,  spoiling  the  amplification.  When  the  voltage  is  removed,  the  Q  is 
restored,  and  oscillation  can  start. 

The  half-wavelength  retardation  system  shown  in  Figure  2-5(b)  is  a  'pulse  on' 
Q-switch.    With  no  voltage  applied  the  Pockels  cell  does  not  alter  the  polariza- 
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Figure  2-6:  Schematic  Cavity  Dump 

tion  of  the  transmitted  beam,  and  the  crossed  polarizers  prevent  oscillation  from 
occurring.  When  switched  on  the  Pockels  cell  shifts  the  polarization  90°  allowing 
the  beam  to  pass  uninterrupted  through  the  two  polarizers. 

The  half-wave  retardation  voltage  for  KDP  is  7500V.  KDP  is  also  water  soluble. 
Normally  the  crystal  is  hermetically  sealed  to  prevent  exposure  to  atmospheric 
moisture,  but  hydrophobic  materials  and  extremely  high  voltages  merit  caution 
when  used  in  a  hydrodynamics  experiment. 


2.5.5.2      Mechanical  Q-Switches 

There  are  two  common  methods  of  mechanically  Q-switching  a  laser.  They  both 
rely  on  mechanically  rotating  a  reflective  element  within  the  laser  cavity.  The  most 
obvious  is  to  rotate  the  mirror  such  that  it  is  only  periodically  parallel  to  the  fixed 
reflector  at  the  opposite  end  of  the  cavity.  This  system  requires  extremely  critical 
alignment  of  the  shaft  upon  which  the  mirror  rotates.  The  alternative  to  rotating 
a  mirror  is  to  use  a  roof  prism,  which  alleviates  the  alignment  requirement. 

Mechanical  switching  is  not  as  flexible  as  electrooptical  switching,  and  is  prob- 
ably not  practicably  capable  of  implementing  the  three  pulse  technique  described 
in  Section  3.1. 
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Figure  2-7:  Cylindrical  Transverse  Electromagnetic  Modes 


2.5.5.3      Cavity  Dumping 

Extremely  short  pulses  can  be  generated  using  a  "cavity  dump".  Cavity  dumps 
allows  the  use  of  two  100%  reflectivity  end  mirrors.  When  the  laser  power  is 
maximized,  the  cavity  dump  redirects  the  beam  outside  of  the  cavity,  rapidly 
dumping  the  entire  optical  energy  of  the  cavity  into  the  output  beam.  The  pulse 
length  is  reduced  to  the  round  trip  transit  time  within  the  cavity.  Cavity  dumps 
can  be  used  to  switch  intermittently  and  continuously  pumped  lasers. 

A  schematic  of  a  cavity  dumped  ruby  laser  is  shown  in  Figure  2-6.  Initially,  no 
voltage  is  applied  to  the  Pockels  cell.  When  the  flash  lamp  is  fired  the  horizontally 
polarized  wave  is  passed  through  the  polarizer  and  not  allowed  to  regenerate. 
Vertically  polarized  light  is  kept  within  the  cavity  and  builds  up  the  vibration  mode 
within  the  cavity.  When  the  energy  stored  in  the  Ruby  has  peaked,  the  Pockels 
cell  is  switched  to  its  half-wave  retardation  voltage.  The  vertically  polarized  light 
is  rotated  to  horizontally  polarized  light  by  the  Pockels  cell  and  passes  through  the 
polarizer  and  out  of  the  cavity.  The  combination  of  the  polarizer,  off- axis  mirror 
and  the  Pockels  cell  essentially  forms  a  voltage  controlled  mirror,  whose  reflectivity 
can  be  change  from  100%  to  0%  and  back  to  100%  during  the  cycle. [6,  8] 
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2.5.6  Transverse  Electromagnetic  Modes 

Laser  cavities  do  not  only  oscillate  longitudinally.  There  are  also  transverse  modes 
of  vibration  which  are  very  similar  to  the  vibration  modes  of  a  drum  head.  Laser 
specifications  typically  include  a  description  of  which  transverse  modes  are  present 
in  the  output.  The  Transverse  Electromagnetic  Modes  (TEM)  are  described  using 
the  symbol  TEM7m?(?.  The  first  two  indices  describe  the  number  of  transverse 
modes,  while  the  third  (often  omitted  in  specification)  describes  the  number  of 
longitudinal  modes.  Figure  2-7  shows  cylindrical  transverse  modes,  along  with  the 
indices  describing  them.  Rectangular  modes  are  also  possible.  It  is  not  difficult 
to  find  TEM oo  ruby  lasers  suitable  for  HPIV  use. 

2.5.7  Output  Frequency  Control 

It  is  possible  to  dramatically  reduce  the  number  of  longitudinal  modes  present  in 
the  output.  The  simplest  method  is  to  use  a  Fabry- Perot  etalon.  This  is  simply 
a  very  small  optical  cavity  that  resonates  at  precisely  the  frequency  desired  in 
the  laser  cavity.  The  etalon  effectively  spoils  the  cavity  Q  for  other  frequencies, 
preventing  them  from  building  up. 

The  simplest  form  of  etalon  is  simply  a  thin  glass  plate.  Internal  reflections 
within  the  plate  selection  interfere  with  the  wave  to  produce  a  sharp  resonance 
peak.  The  effective  thickness  of  the  plate  can  be  changed  by  rotating  the  etalon 
about  an  axis  perpendicular  to  the  cavity  axis. 

Other  etalon  designs  are  possible.  Alternate  designs  include  wedge  shaped 
plates  that  can  be  tuned  by  translating  the  plate,  as  well  as  designs  that  form 
an  optical  cavity  using  an  air  gap  between  two  plates.  The  latter  design  can 
sometimes  be  adjusted  by  varying  the  gas  pressure  or  temperature  within  the 
cavity.  [6,  8] 

2.5.8  Laser  Construction 

Commercial  lasers  for  field  use  are  often  manufactured  in  such  a  way  as  to  limit  the 
ability  of  the  user  to  alter  the  laser.   This  increases  laser  reliability.   Commercial 
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lasers  for  laboratory  use  are  a  different  matter.  As  discussed  above  there  are 
many  additional  components  that  can  be  added  to  the  laser  cavity  to  customize 
and  control  the  laser  output.  Laboratory  lasers  typically  mount  all  components, 
including  the  mirrors  and  laser  rod,  on  a  long  rail  with  additional  space  for  adding 
polarizers,  etalons,  Q-switches  and  cavity  dumps. 
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Chapter  3 


HPIV  System  Design 


3.1      Timing 

The  requirements  to  unambiguously  determine  speeds  less  than  lOm/s  inside  a 
volume  5  cm  on  a  side  impose  severe  restrictions: 

1.  Exposure  times  must  be  short  enough  to  'stop'  the  motion.  In  a  photograph 
the  motion  allowed  during  the  exposure  is  related  to  the  feature  size  of  the 
object.  During  the  photographic  exposure  the  smallest  features  of  an  object 
should  not  move  more  than  approximately  l/10th  their  size  or  the  image  will 
be  unacceptably  blurred.  Likewise  an  acceptable  amount  of  particle  motion 
during  holographic  exposure  has  been  found  to  be  approximately  l/10th  of 
a  particle  diameter  [9]. 

2.  The  interval  between  exposures  must  be  short  enough  to  maintain  the  iden- 
tity of  the  particles,  so  tracking  and  velocity  determination  is  possible.  This 
restricts  the  mean  free  path  length  of  a  particle  between  exposures  to  less 
than  the  distance  between  particles.  There  is  an  additional  problem  of  de- 
termining which  particle  image  is  "oldest',  so  that  direction  of  travel  can  be 
determined. 

3.  The  interval  between  exposures  must  be  short  enough  to  ensure  a  majority 
of  the  particles  captured  in  the  first  exposure  are  within  the  covered  volume 
during  the  following  exposures. 
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For  a  10/iin  particle,  these  requirements  become: 

1.  A  1/is  exposure  results  in  l/10th  diameter  of  travel  during  the  exposure. 

2.  A  10/LiS  exposure  interval  results  in  10  particle  diameters  of  travel  between 
exposures  at  lOm/s. 

The  exposure  interval  also  places  an  upper  limit  on  the  concentration  of  seed 
particles.  Assuring  the  particle  mean  free  path  between  exposures  is  less  than  the 
average  distance  between  particles  requires  that 


Pp 


1  r  particles  fi  particles 

1 %2x  10s- —  =  4  x  10b- = — 

^(lOO^m)3  cm3  in3 


Where  pp  is  the  seed  particle  density,  and  the  mean  free  path  of  the  particles  is 
ten  particle  diameters,  100/im. 

The  direction  problem  is  subtler.  Keeping  the  average  inter-particle  distance 
greater  than  the  mean  free  path  length  between  exposures  significantly  simplifies 
the  task  of  identifying  pairs  of  particle  images  as  being  the  same  particle.  Unfortu- 
nately, even  when  the  two  images  of  a  given  particle  have  been  identified,  there  is  a 
180°  ambiguity  in  the  velocity  of  the  particle.  A  simple  solution  to  this  problem  is 
to  use  three  pulses  rather  than  two.  with  differing  time  intervals  between  pairs  of 
pulses.  For  example,  the  second  l^s  pulse  is  sent  5/is  after  the  first,  then  the  third 
is  sent  10/is  after  the  second.  The  pulse  intervals  are  sufficiently  short  to  permit 
particle  identification,  while  the  difference  in  intervals  allows  determination  of  the 
time  order  of  the  images.  An  additional  advantage  of  the  triple  exposure  tech- 
nique could  be  estimation  of  fluid  acceleration.  This  technique  assumes  the  flow 
speed  will  not  change  appreciably  over  the  distances  traveled  by  particle  during 
the  exposure  sequence.  If  the  particle  were  to  radically  decelerate  to  less  than  half 
of  the  speed  it  traveled  during  the  first  set  of  exposures,  it  would  be  possible  to 
confuse  the  time  order  of  the  exposures. 
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3.2      Recording  and  Digitization 

3.2.1      CCD  Requirements 

The  ultimate  goal  for  computationally  reconstructed  HPIV  is  a  completely  digital 
process  for  determining  the  flow  field  throughout  the  volume.  The  most  likely  can- 
didate for  the  "digital  recording  plane'  is  a  Charge  Coupled  Device  (CCD)  Array. 
Unfortunately,  the  resolution  of  CCD  arrays  is  significantly  less  than  holographic 
film. 

3.2.1.1  Characteristics  of  CCD  arrays 

A  CCD  Array  is  a  semiconductor  architecture  that  transfers  charge  through  stor- 
age areas.  The  CCD  array  is  essentially  an  array  of  charge  "buckets"' .  Photons 
are  captured  by  the  buckets  and  converted  to  charge,  the  buckets  only  perform 
this  conversion  during  a  controllable  integration  period.  More  photons,  or  higher 
energy  photons,  hitting  the  array  during  the  integration  period  result  in  a  higher 
element  charge.  After  integration  the  CCD  Array  support  electronics  scan  through 
each  bucket  and  report  the  charge  level  in  each  bucket.  This  charge  is  proportional 
to  the  intensity  of  the  image  focused  on  the  array.  [10] 

For  holographic  applications  the  fundamental  characteristic  of  CCD  Arrays 
that  differentiates  them  from  film  is  the  limited  number  of  recording  elements. 
High  quality  photographic  film  can  have  grains  smaller  than  1/mi  over  an  area  12- 
20  centimeters  square,  or  better  than  120,000x120,000  recording  elements.  The 
best  CCD  arrays  currently  commercially  available  have  element  sizes  on  the  order 
of  10/mi,  but  the  total  array  size  is  only  about  5  centimeters  on  a  side,  yielding 
an  array  of  approximately  5000x5000  elements.  This  is  the  root  of  the  limitation 
in  holographic  applications. 

3.2.1.2  CCD  Resolution  Requirements 

3.2.1.2.1  Diffraction  patterns  of  circular  objects  Parrent  and  Thomp- 
son analytically  derived  the  diffraction  pattern  for  an  opaque  circular  object  in  a 
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coherent  field.  [11]  For  a  circular  particle,  the  field  and  intensity,  respectively,  are: 

ka2         (  kr2\        ( k.ar 

ib(r)     =    coskz sin     kz  +  ——     Ai 

r  \  Iz 
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Where,  r  is  measured  from  the  center  of  the  object  projected  into  the  recording 
plane,  and  A(q)  =  Ji(a)/a,  a  Bessel  function  of  the  first  kind  divided  by  its 
argument. 

An  illustrative  example  is  plotted  in  Figures  3-1  and  3-2.  The  parameters 
more  characteristic  of  HPIV  systems  result  in  plots  that  do  not  acceptably  show 
the  salient  features  of  the  intensity  distribution. 

3.2.1.2.2  Recording  the  Diffraction  Pattern  Equation  (3.1)  describes  the 
diffraction  pattern  recorded  for  a  single  particle  with  circular  cross  section.  All 
previous  discussion  has  assumed  this  diffraction  pattern  is  perfectly  recorded.  In 
reality  the  recording  medium  has  limited  bandwidth  (i.  e.  a  grain),  and  limited 
physical  extent. 

Figure  3- 1  shows  a  high  frequency  sine  function  modulated  by  a  lower  frequency 
A(q')  envelope  function.  The  envelope  function  is  dependent  on  both  particle 
diameter  and  particle  distance  from  the  recording  plane.  The  high  frequency 
sine  function  is  dependent  only  on  particle  distance.  Accurately  recording  the 
high  frequency  component  is  the  most  critical  factor  in  obtaining  good  depth 
discrimination. 

Vikram[3]  recommends  recording  the  first  three  side  lobes  formed  by  the  bessel 
envelope  function.    The  spacing  of  the  fine  fringes,  Ax,  in  the  third  lobe  can  be 
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Figure  3-2:  Enlarged  plot  of  the  third  side  lobe  in  Figure  3-1 


found  using  the  argument  of  the  sine  term  in  (3.1),  kr2/2z: 


U    2 
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~2z~ 


=     2tt 
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2  1 
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(3-2) 


Equation  (3.2)  should  be  evaluated  in  the  region  containing  the  third  side  lobe. 
The  mid-radius  of  the  third  side  lobe,  f3  is: 


Ji{an)/an    =    0     (n  =  1,2,  ...,oo) 


o, 


r3 


an  =  3.832,  7.06,  10.173,  13.323,  16.470. 
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(3-3) 


Equations  (3.3)  and  (3.2)  give  f3  =  0.189cm  and  Ax  =  3.65  x  10-3cm  for 
the  parameters  plotted  in  Figures  3-1  and  3-2.  Using  z  =  5cm  and  a  =  10/rni, 
with  A  =  690nm.  more  reasonable  for  HPIV  purposes,  f3  =  0.472cm  and  Ax  = 
7.31  x  10-4cm. 

The  parameter  f3  dictates  the  minimum  size  of  the  hologram,  with  no  mag- 
nification, and  is  easily  met.    The  parameter  Ax  dictates  the  minimum  grain  of 
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the  hologram  recording  medium.  Since  enlargement  or  reduction  of  the  hologram 
during  recording  is  easily  accomplished  using  lenses,  the  absolute  size  of  the  fringe 
spacing  is  not  important.  Rather,  the  important  quantity  is  the  ratio  the  fringe 
spacing  to  the  total  size  of  the  system,  so  that  the  fine  fringes  produced  by  par- 
ticles at  opposite  extremes  of  the  volume  can  be  resolved.  The  sampling  theorem 
suggests  that  the  recording  medium  should  be  able  to  sample  at  twice  the  high- 
est spatial  frequency  required,  or  1/2  the  fringe  spacing.  To  achieve  the  desired 
accuracy,  a  CCD  array  would  have  to  have 

2  •  0.05cm  nr.nn 

13680  pixels  square. 


7.31  x  10-4cm 


This  resolution  will  adequately  record  the  particle  diffraction  patterns.  Unfortu- 
nately no  CCD  arrays  are  commercially  available  that  reach  this  resolution.  The 
best  currently  available  are  in  the  range  of  5000  pixels  square. 

3.2.2      Exposure  Requirements 

While  the  performance  characteristics  of  CCD  arrays  are  well  known  it  is  difficult 
to  predict  the  amount  of  power  incident  on  the  array  from  the  suspended  particles. 

The  resolution  analysis  was  performed  assuming  the  finest  necessary  fringes 
were  detected.  The  incident  power  depends  heavily  on  the  power  output  of  the 
laser,  the  absorption  of  the  fluid,  and  the  reflectivity  of  the  particles,  and  the 
viewing  window  material. 

The  exposure  problem  is  different  in  character  from  the  resolution  problem. 
The  resolution  problem  has  to  wait  for  better  CCD  technology,  or  lower  the  ac- 
curacy requirements  for  HPIV  significantly.  The  exposure  problem  can  be  solved 
many  ways,  the  simplest  (but  not  necessarily  the  least  expensive)  way  is  to  increase 
the  laser  power. 


36 


Chapter  4 


Simulation 


4.1  Motivation 

A  primary  factor  preventing  HPIV  from  becoming  an  easily  and  widely  used  tech- 
nique is  the  data  reduction  effort.  As  previously  mentioned,  physical  reconstruc- 
tion techniques  are  simple,  but  require  costly  equipment,  and  inordinate  amounts 
of  time.  If  computational  reconstruction  methods  can  be  perfected,  HPIV  could 
be  become  more  practical. 

There  are  several  possibilities  for  computational  reconstruction  algorithms, 
one  of  which  (Onural's  method)  is  presented  in  Section  5.1.  The  methods  must 
be  evaluated  with  realistic  input  in  order  to  gauge  their  effectiveness.  Obtaining 
realistic  input  is  difficult.  There  have  been  experiments  done  using  holographic 
techniques  to  study  small  ocean  biological  specimens,  but  obtaining  the  actual 
holograms  produced  is  difficult.  An  additional  difficulty  is  that  you  don't  know 
precisely  what  was  contained  in  the  original  image  so  calibrating  the  reconstruction 
algorithm  is  nearly  impossible. 

4.2  Hologram  Synthesis 

There  are  several  methods  applicable  to  synthesizing  holograms  similar  to  those 
expected  from  HPIV.  The  most  straightforward  methods  calculate  the  individual 
contribution  from  each  particle  to  each  pixel  on  the  simulated  recording  plane. 
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These  methods  are  the  most  accurate  and  flexible  but  suffer  from  excessive  com- 
putational requirements  when  large  numbers  of  particles  are  being  simulated. 
One  way  around  this  problem  is  hinted  at  by  Equation  (2.18): 


ii>{x,  y)  =  ^2  Un(x'  y)  ®  hz»  (x'  2/) 


We  can  subdivide  the  region  being  simulated  into  a  large  number  of  thin  slabs. 
Within  each  slab  there  can  be  a  large  number  of  particles,  but  each  will  be  treated 
as  if  it  had  the  same  z  coordinate.  We  then  calculate  the  contribution  of  each 
slab  to  the  total  hologram  rather  than  the  contribution  from  each  particle.  With 
a  large  number  of  particles  we  can  achieve  significant  performance  increases  with 
little  degradation  of  quality. 

Evaluating  (2.18)  is  quite  efficient.  Using  the  convolution  theorem: 


V;  {x,  y)  =  ^2  u»  (Xi  y)  ®  h**  (^  y) 
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Where UZn  {ujT,uy)  =  T[U  (x,ij))  and  UZn  {ux,uy)  =  F[hZri  (x,y)].  ButHZn  (ux,uy) 
can  be  analytically  evaluated  [12], 


UZn  {ux,ujy)  =  T[hZn]  =  exp 


•      Zn    (     2    i        2\ 
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yielding. 


ip{x,y)  =T  l 


/Mn  (^T^y)  ■  exp 


■  ^zn    (     2    i        2\ 


(4.3) 


The  algorithm  to  evaluate  (4.3)  is  shown  schematically  in  4-1.  and  computer 
code  to  perform  the  algorithm  is  given  in  Sections  A.  1.1  and  A.  1.2. 
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Figure  4-1:  Hologram  synthesis  algorithm 
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4.3     HoloSynth  Simulator 
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Figure  4-2:  HoloSynth  user  interface 


4.3.1      Introduction 


HoloSynth  assists  in  synthesizing  holograms  similar  to  those  obtained  from  Holo- 
graphic Particle  Image  Velocimetry.  It  provides  the  user  with  a  simple  way  to 
place  any  number  of  particles  into  a  cubic  volume  in  preparation  for  generating 
a  synthetic  hologram.  The  user  can  manually  place  particles  anywhere  in  the  3- 
dimensional  region  as  well  as  fill  arbitrary  rectangular  solid  regions  with  randomly 
distributed  particles. 

HoloSynth  uses  the  method  for  synthesizing  holograms  described  in  Section 
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Figure  4-3:  HoloSynth  main  controls 

4.2.  In  addition  it  prepares  an  input  file  for  an  alternate  hologram  synthesizer 
so  that  different  synthesis  algorithms  may  be  easily  compared,  without  modifying 
the  HoloSynth  source  code.  The  output  file  is  also  useful  in  analyzing  the  results 
of  reconstruction  algorithms. 

HoloSynth  was  developed  under  Red  Hat(R)Linux  version  5.1  and  tested  under 
Microsoft  Windows  NT®. 

4.3.2      The  Interface 

4.3.2.1      The  Main  Controls 

Figure  4-3  shows  the  main  controls  for  HoloSynth. 

Each  button  will  open  up  a  new  window  that  controls  a  particular  aspect  of 
HoloSynth.  The  windows  are: 

1.  3-D  View  of  the  current  particle  configuration. 

2.  Plane  View,  which  allows  manual  placement  of  particles  on  an  arbitrary 
plane. 

3.  Select  Region,  which  allows  selection  of  an  arbitrary  rectangular  solid  region 
within  the  view  space.  From  this  control  panel  you  can: 


•  Clear  all  or  any  part  of  the  volume. 
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□    HoloSynth  1.0: 3-D  View 
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Figure  4-4:  3-D  View  window 

•  Fill  all  or  any  part  of  the  volume  with  randomly  placed  particles. 

4.  The  Synthesizer  Settings  window  which  controls  the  physical  characteristics 
of  the  simulation  used  for  the  hologram  synthesis. 

Through  the  file  menu  the  user  can  open  and  save  particle  fields,  and  exit 
HoloSynth. 


4.3.2.2      The  3-D  View 

The  3-D  View,  shown  in  Figure  4-4,  provides  a  three  dimensional  view  of  the 
current  particle  configuration.  There  are  three  basic  functions  controlled  by  this 
window:  the  magnification  (zoom)  of  the  view,  the  angle  and  translation  of  the 
view,  and  the  position  of  the  active  and  two  clipping  planes.  No  direct  editing 
takes  place  in  this  window,  it  is  merely  a  visualization  aid. 

The  viewing  "camera  position"  is  controlled  by  the  rollers  and  sliders  to  the  left 
and  below  the  viewing  window.  The  slider  above  the  viewing  window  controlling 
the  magnification  of  the  camera. 

The  various  view  control  planes  are  controlled  using  the  slide  bars  to  the  right  of 
the  3-D  view  window.  They  appear  as  colored  squares  drawn  around  the  periphery 
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a    HoloSynth  1.0:  Active  Plane  View 
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Figure  4-5:  HoloSynth  plane  editor 

of  the  cube.  Only  the  active  plane  is  visible  in  Figure  4-4,  the  clipping  planes  are 
hidden  at  the  ends  of  their  travel.  On  a  color  display  the  active  plane  is  bright 
green,  the  rear  cliiping  plane  is  blue,  and  the  front  clipping  plane  is  magenta.  The 
active  plane  determines  the  z-coordinate  that  the  plane  editor  uses  when  adding 
new  particles.  The  front  and  rear  clipping  plane  restrict  the  range  of  z  coordinates 
visible  in  both  the  plane  view  and  the  3-D  view.  In  all  cases,  particles  lying  on 
the  active  plane  are  visible  within  both  views.  If  the  user  wants  to  view  only  the 
particles  on  the  active  plane,  then  he  places  both  clipping  planes  at  the  same  limit 
of  travel. 


4.3.2.3     The  Plane  View 

The  plane  view  projects  all  particles  visible  in  the  3-D  view  onto  a  single  plane. 
Particles  that  are  on  the  active  (green)  plane  in  the  3-D  view  are  highlighted  green 
in  the  plane  view.  Particles  can  be  manually  positioned  or  deleted  on  the  active 
plane. 

To  add  a  particle,  simply  left  click  within  the  black  portion  of  the  active  plane 
view.  Accurate  positioning  of  particles  can  be  performed  in  two  ways: 
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HoloSynth  1.0: 3-D  Selection 
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Figure  4-6:  3-D  region  selection 


1.  Turn  on  the  coordinate  grid.  At  the  bottom  of  the  plane  view  window  is  a 
selection  list  allowing  presentation  of  grids  of  various  densitites.  High  grid 
densities  are  convenient  at  high  magnifications.  While  a  grid  is  displayed 
the  particles  automatically  snap  to  the  nearest  grid  intersection.  This  can 
be  controlled  by  the  Display->Snap  to  Grid  menu  item,  and  can  also  be 
toggled  using  the  'g'  key. 

2.  Use  the  real-time  coordinate  diplay  Any  time  the  mouse  is  over  the  edit 
window  its  complete  3-D  coordinate  is  show  at  the  bottom  of  the  window. 

Deletion  is  accomplished  by  dragging  with  the  right  mouse  button  over  the 
region  containing  the  particles  to  be  deleted.  A  red  selection  border  will  appear 
showing  the  selected  region.  To  delete  inside  the  selected  region,  use  the  menu 
item.  Edit->Delete  Inside  Selection  or  the  delete  key.  Using  the  Edit->Delete 
Outside  Selection  menu  item  will  clear  all  particles  visible  in  the  plane  editor 
but  outside  the  selection  box.  The  region  deleted  is  bounded  in  the  x  and  y  di- 
rections by  the  selection  box,  and  in  the  z  dimension  by  the  clipping  planes.  If 
a  particle  is  visible  in  the  plane  editor,  regardless  of  whether  it  is  on  the  active 
plane,  it  can  be  deleted. 
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4.3.2.4  Select  Regions 

The  3-D  Selection  window,  shown  in  Figure  4-6,  allows  the  user  to  define  a  solid 
rectangular  sub- volume  within  the  3-D  view  and  perform  the  following  operations: 

1 .  Clear  the  entire  volume. 

2.  Clear  all  particles  inside  the  selected  sub- volume 

3.  Clear  all  particles  outside  the  selected  sub- volume 

4.  Fill  the  entire  volume  with  randomly  placed  particles.  The  slider  at  the  right 
of  the  window  controls  how  many  particles  will  be  added. 

5.  Fill  the  selected  sub- volume  with  randomly  placed  particles.   The  slider  at 
the  right  of  the  window  controls  how  many  particles  will  be  added. 

6.  Fill  outside  the  selected  sub-volume  with  randomly  placed  particles.    The 
slider  at  the  right  of  the  window  controls  how  many  particles  will  be  added. 

The  sub-volume  is  displayed  as  a  translucent  yellow  box  with  the  3-D  view.  Each 
wall  of  the  selection  region  is  controlled  by  its  own  slider.  The  region  will  only  be 
displayed  in  the  3-D  view  if  the  Show  Selection  Box  button  has  been  selected. 

4.3.2.5  Hologram  Synthesis 

Figure  4-7  shows  the  synthesizer  controls.  This  allows  easy  control  of  the  physical 
size  of  the  simulation  as  well  as  recording  plane  distance  and  number  of  planes  used 
to  subdivide  the  simulation  volume.  The  generated  hologram  will  be  displayed  in  a 
new  window  when  the  computations  have  completed.  If  the  product  is  satisfactory, 
the  user  may  save  the  hologram  to  a  file  using  the  Truevision  Targa  format  or  the 
Adobe  Photshop  Raw  format. 

4.3.3      Requirements 

HoloSynth  is  written  completely  in  C++.  The  Graphical  User  Interface  uses 
the  freely  available  platform  independent  user  interface  toolkit  FLTK[13].    The 
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Figure  4-7:  Synthesizer  settings 

hologram  synthesizer  use  the  FFTW  fast  fourier  transform  library  available  freely 
from  MIT [14].  The  program  runs  under  X- Windows,  and  Microsoft  Windows  with 
OpenGL[15]  or  the  Mesa  3-D  graphics  library  (an  OpenGL  functional  equivalent)  [16]. 


4.3.4     Program  Architecture 

At  its  heart  HoloSynth  is  a  fancy  list  editor,  editing  a  list  of  particle  coordinates 
and  radii.  The  core  data  structure  is,  appropriately,  a  linked  list  of  particle  objects. 
This  allows  an  arbitrary  number  of  particles  to  be  added,  restricted  only  by  the 
memory  available. 

Figure  4-8  shows  the  schematic  relationship  between  all  program  modules  that 
make  up  HoloSynth.  The  complete  source  code  to  HoloSynth  is  presented  in 
Appendix  A.l. 
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Figure  4-8:  HoloSynth  architecture 
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Chapter  5 


Computational  Reconstruction 


5.1      Onural's  Method 

In  1985,  Levent  Onural  proposed  a  technique  based  on  a  truncated  inverse  filter [12, 
17].    The  method  is  summarized  here,  source  code  is  presented  to  perform  the 
algorithm  in  Sections  A. 2.1  and  A. 2. 2. 
Starting  from  (2.17), 

ipz(x,  y)  =  U(x,  y)  <8>  hz(x,  y) 

We  let  U  (x,y)  =  1  —  a  (x,  y)  where  a  (x,  y)  is  the  opacity  of  the  given  plane.  Now 
write  the  recording  plane  intensity, 

h  (x,  y)    =    ipz  (x,  y)  'ip*z  (x,  y) 

=     \{l-a{x,y))®hz(x,y)\2 

=     1  -  a   (.t,  y)  ®  h*z  (x,  y)  -  a  (x,  y)  ®  hz  ( x,  y) 

+  |a(x,2/)(g)/iz(x,2/)|2  (5.1) 

=     1  -  a  (x,  ?/)  <8>  [/i*  (x,  y)  +  hz  (x,  t/)] 
=    l-a(x,j/)®2Re{/iI(x,2/)}  (5.2) 

Where  |a(x,?/)  0  hz  (x,y)\    was  dropped  as  insignificant  in  (5.1).   Equation  (5.2) 
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shows  the  intensity  modeled  as  a  DC  shifted  linear  system  with  impulse  response. 

g{x,y)  =  2Re{hz{x,y)}  (5.3) 

No  inverse  of  g  (x,  y)  in  the  normal  sense  exists,  However,  Onural  showed  in 
[12]  that  a  series  approximation  to  the  Fourier  transform  of  g  can  be  found.  If 
G  {ujx,ujy)  =  T  [g  (x,  y)}  then, 

Ci,K.^)-c«^(^  +  ^)+    J^       tf-'M|l08i'-1     S,(M,>,)     (5.4) 

71  fc=l  ^  ' 

where  pog2  A:J  is  the  integer  obtained  by  truncating  the  value  of  log2  k  and, 

Sk  (ux,uy)  =  (-l)k cos  ^  (2*  +  1)  (u£  +  wj)  (5.5) 

yielding. 

^M1(^2/)=^""1[^M(^^y)]  (5-6) 

Looking  back  at  the  hologram  recording  given  by  (5.2), 

a(x,y)®2Re{hz{x,y)}    =    l-I{x,y) 
T[a{x,y)]Q(ux,ujy)    =    T  [1  -  /  (x,  y)] 

az{x,y)    =    T-l[T[l-I{x,y))g^z{ux,uy)\        (5.7) 

yielding  the  original  opacity  function. 

5.2      Simulation  Results 

5.2.1      Synthesis  and  Computational  Reconstruction 

Figure  5-2  is  a  portion  of  a  hologram  generated  by  HoloSynth.  The  particle  config- 
uration consisted  of  2328  particles.  2000  of  the  particles  were  randomly  distributed 


5.2.     SIMULATION  RESULTS 


49 


Figure  5-1:  Schematic  of  the  particle  configuration  used  to  generate  Figure  5-2 


throughout  a  5cm  cube.  Along  the  center  plane  parallel  to  the  wavefront  328  par- 
ticles were  positioned  to  spell  the  letters  "MIT".  Figure  5-1  illustrates  a  schematic 
of  that  configuration.  The  synthetic  hologram  itself  was  2048x2048  pixels,  and 
the  volume  was  subdivided  into  1024  planes,  for  a  z-axis  discretization  of  0.04mm. 
The  simulated  particle  diameter  was  25/um.  The  portion  of  the  hologram  shown  in 
5-2  is  512x512  pixels  and  its  upper  left  corner  is  205  pixels  from  the  left  extremity 
of  the  original  and  741  pixels  from  the  top.  Although  it  contains  the  'M'  from 
"MIT"  .  it  is  impossible  to  distinguish  from  the  diffraction  patterns  produced  by 
the  other  particles. 

Figure  5-3  through  5-8  show  the  results  of  reconstructing  Figure  5-2  at  various 
distances  between  30.0cm  and  30.5cm  from  the  recording  plane  using  Onural's 
method  with  M  =  5.  The  images  generated  are  the  opacity  functions,  with  black 
being  transparent,  white  being  opaque. 

In  Figure  5-3(a)  the  'M1  is  clearly  visible  on  a  computer  monitor  with  some 
slight  clouds  from  out-of-focus  particles.  Figure  5-4  is  the  same  hologram  recon- 
structed at  30.1cm.  Here  there  is  little  visible  difference  from  the  previous  figure. 
This  is  an  artifact  of  the  way  the  decoder  transform  real  valued  intensities  into 
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gray  scale  bitmap  images. 

To  translate  from  an  arbitrary  field  of  real  number  to  a  field  of  8  bit  integers, 
the  decoder  first  finds  the  maximum  and  minimum  values  present  in  the  real  array, 
it  then  scales  and  translates  each  of  these  values  to  fill  the  range  from  0  to  256, 
in  order  to  maximize  use  of  dynamic  range  available  in  an  8  bit  gray  scale  image. 

The  30.1cm  reconstruction  has  a  much  lower  maximum  intensity  than  the 
30.0  cm  reconstruction.  Volume  reconstruction,  where  many  planes  will  be  recon- 
structed at  a  time,  in  order  to  reconstruct  flow  fields,  will  have  to  keep  track  of 
this  scaling  issue  between  planes. 

The  important  thing  to  note  about  this  series  of  decodes,  is  that  while  the 
source  hologram  was  constructed  using  a  z  discretization  of  0.04mm,  the  recon- 
struction algorithm  2-discrimination  was  of  the  order  of  millimeters.  Thus,  the 
original  image  could  have  been  discretized  into  64  steps,  or  0.78mm  slabs,  requiring 
l/16th  the  computation  time  with  little  or  no  loss  of  reconstruction  accuracy. 

5.2.2      Image  Postprocessing 

The  next  step  in  reconstruction  is  to  automatically  detect  particles  in  each  plane 
and  register  their  x,  y,  and  z  coordinates,  rather  than  simply  displaying  an  im- 
age. Repeatedly  reconstructing  at  different  depths  and  detecting  the  particles  will 
replace  the  time  consuming  3-D  scanning. 

In  Figures  5-3  through  5-8  subfigure  (b)  and  (c)  show  the  results  of  applying 
some  simple  image  processing  techniques  to  the  original  reconstructions  to  enhance 
the  output.  Subfigure  (b)  in  each  figure  applies  'edge  detection',  which  simply 
filters  the  image  marking  any  region  where  the  gradient  of  the  image  exceeds 
a  certain  threshhold.  In  the  case  of  the  grayscale  images  being  used  here,  the 
gradient  is  simply  the  difference  in  grayscale  value  between  adjacent  pixels.  [18] 
Subfigure  (c)  in  each  figure  'threshholds'  the  image,  changing  a  pixel  to  white 
(value  255)  or  black  (value  0)  depending  on  whether  or  not  the  pixel  value  exceeds 
a  preset  threshhold.1 


'The  image  processing  was  accomplished  with  the  freely  available  GNU  Image  Manipulation 
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The  combination  of  edge  detection  and  threshholding  can  be  used  very  effec- 
tively to  automatically  detect  particle  positions.  At  the  end  of  the  image  process- 
ing, there  are  only  regions  of  black  (value  0)  with  very  small  regions  of  white  (value 
255).  The  white  regions  correspond  to  particles  and  can  easily  be  automatically 
tabulated  and  mapped  into  x  and  y  coordinates. 


Program  (GIMP,  http://www.gimp.org),  which  is  very  similar  to  Adobe  Photoshop®. 
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Figure  5-2:  Portion  of  a  synthetic  hologram.  The  original  synthetic  hologram  was 
2048x2048  pixels.  This  portion  is  512x512  pixels,  cropped  from  the  orig- 
inal hologram  to  make  the  images  presentable  in  printed  form. 
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(a)  Raw  Reconstructed  Image 


(b)  Image  after  enhancing  with  edge  de- 


tection. 


(c)  Image  after  enhancing  with  edge  de- 
tection and  threshholding. 


Figure  5-3:  Figure  5-2  decoded  at  a  depth  of  0.300m. 
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(a)  Raw  Reconstructed  Image 


(b)  Image  after  enhancing  with  edge  de- 
tection. 


(c)  Image  after  enhancing  with  edge  de- 
tection and  threshholding. 


Figure  5-4:  Figure  5-2  decoded  at  a  depth  of  0.301m. 
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(a)  Raw  Reconstructed  Image 


(b)  Image  after  enhancing  with  edge  de- 
tection. 


(c)  Image  after  enhancing  with  edge  de- 
tection and  threshholding. 


Figure  5-5:  Figure  5-2  decoded  at  a  depth  of  0.302m. 
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(a)  Raw  Reconstructed  Image 


(b)  Image  after  enhancing  with  edge  de- 
tection. 


(c)  Image  after  enhancing  with  edge  de- 
tection and  threshholding. 


Figure  5-6:  Figure  5-2  decoded  at  a  depth  of  0.303m. 
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(a)  Raw  Reconstructed  Image 


(b)  Image  after  enhancing  with  edge  de- 
tection. 


(c)  Image  after  enhancing  with  edge  de- 
tection and  threshholding. 


Figure  5-7:  Figure  5-2  decoded  at,  a  depth  of  0.304m. 
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(a)  Raw  Reconstructed  Image 


(b)  Image  after  enhancing  with  edge  de- 
tection. 


(c)  Image  after  enhancing  with  edge  de- 
tection and  threshholding. 


Figure  5-8:  Figure  5-2  decoded  at,  a  depth  of  0.305m. 
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Chapter  6 

Conclusions  and 
Recommendations 


It  is  clear  that  computational  reconstruction  of  HPIV  is  practicable.  The  proce- 
dure that  Onural  proposed  in  1985  was  not  then  practical  due  to  the  computational 
power  commonly  available.  Modern  computing  can  easily  handle  realistic  holo- 
grams and  the  numeric  techniques  of  normal  PIV  could  be  extended  to  deal  with 
the  particle  identification  problem  in  HPIV.  The  simulator  presented  in  this  thesis 
will  be  useful  in  any  setting  when  attempting  to  validate  particular  reconstruction 
algorithms. 

Future  work  should  concentrate  on  extending  the  reconstruction  prosses  fully 
into  three  dimensions  and  developing  particle  identification  techniques.  Automat- 
ically determining  the  velocity  vector  for  a  particle  is  the  ultimate  goal. 
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Appendix  A 


Source  code 


A.l      HoloSynth:  Hologram  Synthesizer. 

This  is  the  source  code  for  HoloSynth,  a  program  for  easily  creating  particle  configuration 
and  holograms  for  use  calibrating  computational  reconstruction  algorithms. 

A. 1.1      OnuralSynthesizer.hpp 

#ifndef  ONURALSYNTHESIZER_H 

//  .'include  guard 

#def  ine  ONURALSYNTHESIZERJi  1 

#include  "Synthesizer  .ripp'1 

#include  "ParticleList  .hpp"' 

#ifndef  FLMATH 

#def  ine  FLMATH 

#include  <FL/math.h> 

#endif 

#include  "RGBImage .hpp" 

/*f\class  OnuralSynthesizer 
* 

*\brief  The  algorithm  use  for  actually  generating  the  hologram, 
* 

*  The  forward  diffraction  algorithm  implement  a  method  described  in 

*  "Digital  Decoding  of  In- Line  Holograms''  a  PhD  thesis  by  Leven  Onural 

*  */ 

class  OnuralSvnthesizenpublic  Synthesizer 

{ 
public: 

///  created  with  a  reference  to  the  current  particle  list 
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OnuralSynthesizer  (PartideList  *pl); 

~OnuralSynthesizer(); 

///  Generate,  organized  all  the  FFTs  and  diffraction  calculations. 

int  generate(); 

private: 

RGBImage  *hologram; 

void  diffract, (fftw .complex  *image); 
void  clearField(fftw_complex  *field,  fftw_real  val); 
void  plot,Particle(particle  *c.fftw_complex  *plane); 
void  plot  (int  p_X,  int  p_Y,  fftw_complex  *plane); 
void  addto(fftw .complex  *al,  fftw .complex  *a2); 
void  MagtoR.GB(fftw_complex  *data); 

void  cleanup(); 

fftwnd_plan  fftwFwdPlan; 
fftwnd_plan  fftwBwdPlan; 


fftw_complex  *oplane; 
fftw _complex  *hplane; 


}; 


#endif 

A. 1.2      OnuralSynthesizer.hpp 

#include  <FL/F1.H> 
#include  "ProgressUI  .hpp'1 
#include  " Scroll ImageUI  .hpp" 
#include  <FL/Fl_Image  .H> 
#include  "  OnuralSynthesizer  .  hpp'" 

#include  ''StatusBox.hpp" 

#define  MIN(a,  b)  (((a)  <  (b))  ?  (a)  :  (b)) 

#def  ine  MAX(a.  b)  (((a)  >  (b))  ?  (a)  :   (b)) 


OnuralSynthesizer: :OnuralSynthesizer(ParticleList  *pl) 

:Synthesizer(pl) 
{ 
! 

OnuralSynthesizer:  :~OnuralSynthesizer() 

{ 

cleanup(); 
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\ 


int  OnuralSynthesizer::generate() 

{ 

int  progresscount=0: 

progressUI— >setMax(5+5*d_Planes); 
progressUI— mpdate(O); 
progressUI— Preset  ( ) : 
progressUI  — >show( ) ; 
Fl::wait(0.0); 

particle  *curr; 

int  pixels,  il,  pcount,  Ncalc,  Mcalc; 

fftw.real  zgrain,  Zmin,  Zmax; 


*  Initialization. 

pixels=2*(2*N/2+l)*2*M; 

Ncalc=2*N;  // Double  the  size  for  padding  out.  the  edge  effects. 
Mcalc=2*M; 

if(progressUI— >set,Status( "Allocating  memory.  .  .")){ 
cleanup(); 
progressUI— >hide(); 
return  0; 

}; 

hplane=(fftw .complex  *)  calloc(  pixels,sizeof(fft,w_complex)); 
oplane=(fftw_complex  *)  malloc(  pixels*sizeof(fftw_complex)); 

progressUI— >update(++progresscount); 

zgrain=(rTtw_real)l/(2*d_Planes);//77ie  thickness  of  a  HALF  plane 

if(progressUI->setStatus("Generating  FFTW  Plans.  .. ")){ 
cleanup(); 
progressUI— »hide(); 
return  0; 

}; 


fftwFwdPlan=fftw2d_create_plan(Ncalc,  Mcalc,  FFTW_FORWARD. 
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FFTWJN.PL  ACE); 

fftwBwdPlan=fftw2d_create_plan(Ncalc,  Mcalc,  FFTW_BACKWARD. 
FFTW_IN_PLACE); 

progressUI^update(++progresscount); 

// Remember  we  will  eventually  need  to  split  out  the  vector  end  of  things. 

pList— ►sortListOnZQ; 


curr=pList— >first(); 
for(il=l;il<=d_Planes;il++){ 

Fl::wait(0); 

if(progressUI— >setStatus("  Clearing  Object  Plane... ")){ 
cleanup(); 
progressUI^hide( ) ; 
return  0; 

}; 

clearField(oplane,  1.0); 

progressUI— >update(++progresscount); 

Zmin=(il-l)*2*zgrain-0.5; 
Z=(Zmin+zgrain)*depth;//P/?.yszca/  Z 
Zmax=Zmin+2*zgrain; 

pcount=0; 

if(progressUI— >-setStatus(  "Plotting  Particles.  .  .")){ 
cleanup(); 
progressUI— »hide(); 
return  0: 

}; 

while((curr!=NULL)&&(curr->Z>Zmin)&&(curr->Z<=Zmax)){ 
plotParticle(curr.oplane) ; 
pcount++; 
curr=curr— >next; 

} 

progressUI— »update(++progresscount); 

if(pcount>0){ 

if(l)rogressUI->setStatus("Performing  Forward  FFT...")){ 
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cleanup(); 
progressUI— diide(); 
return  0; 

}; 

fftwnd_one(ff'twFwdPlan,  oplane,  0); 
progressUI— >update(++progresscount); 

if(progressUI— >setStatus( "Diffracting.  .  . ")){ 
cleanup(); 
progressUI— >-hide( ) ; 
return  0; 

}; 

diffract  (oplane); 

progressUI— >update(++progresscount); 

if(progressUI— »setStatus( "Adding.  .  . ")){ 
cleanup  (); 
progressUI— diide(); 
return  0; 

}; 

addto(hplane,  oplane) ; 

progressUI— mpdate(++progresscount); 

} 

progresscount=5*il+2; 

progressUI— mpdate(progresscount); 

} 

hplane[0].re=  1000000; 

progressUI— »setStatus("  Performing  Backward  FFT..."); 
progressUI— >update(++progresscount); 
fftwnd_one(fftwBwdPlan,  hplane,  0); 
progressUI^setStatus("Calculating  Intensity.  .  ."); 

MagtoRGB  (hplane); 

J  J  The  hologram  is  finished,  now  display  it... 

FUmage  *image=new  FUmage( hologram— »data(),  N,  M); 

progressUI->set,Status(  "Displaying.  .  ."); 

progressUI— mpdate(++progresscount); 

ScrollImageUI  *si=new  ScrollImageUI( image  ,  0,0,  500,500); 

si— >show(); 

progressUI^diideQ; 
cleanup(); 

return  1: 
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void  OnuralSynthesizer::addto(fftw .complex  *a,  fftw .complex  *b) 

{ 

int  il; 

for  (il=0;  il  <  (  4  *  N  *  M  );  il++){ 

if((fabs(b[il].re)>=0.0000001)&&(fabs(b[il].im)>=0.0000001)) 

{ 

a[il].re  +  =  b[il].re: 
a[il].im  +=  b[il].im; 

} 
} 
I 

void  OnuralSynthesizer::clearField(fftw_complex  *field,  fftw_real  val) 
{ 

int  il; 

for(il=0;  il  <  (  4  *  N  *  M  );  il++){ 

field  [il]. re  =  val; 

field[il].im  =  0.0; 
} 

} 

void  OnuralSynthesizer::plotParticle(particle  *c,fft,w_complex  *plane) 

{ 

int  x,  y,  r; 

x=(int)((c-»X  +  0.5)*N); 

y=(int)((0.5  -  c->Y)*M); 
r=(int)(c-»R*N); 

switch  (r) 

{ 

case  6: 

plot(x-l,  y+1,  plane); 

plot(x-l,  y-1,  plane); 

plot(x+l,  y-1,  plane); 

plot(x+l,  y+1.  plane); 
case  2: 
case  3: 
case  4: 
case  5: 

plot(x,  y-1,  plane); 

plot(x,  y+1,  plane); 

plot(x-l,  y,  plane); 

plot(x+l,  y,  plane); 
default: 

plot(x,  y,  plane); 
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} 

void  OnuralSynthesizer::plot(int  p_X,  int  p_Y,  fftw_complex  *plane) 

{ 

int  pos=2  *  p_Y  *  N  4-  p_X; 

if(pos>0){ 

plane[  pos].re=0.0; 

plane[  pos].im=0.0; 

} 

} 

#def  ine  F_PI  4.0*3.1415926535 

void  OnuralSynthesizer::diffract(fftw_complex  *image) 

{ 

int  x.  y,  Ncalc,  Mcalc; 
long  ik; 
fftw_real  lz.  cc: 
fftw_real  wx,  wy,  a,  ca,  sa; 

Ncalc=2*N; 
Mcalc=2*M; 

lz=L*(distance-Z); 
ik=0: 


/*  FFTW  puts  frequency  space  origin  at  the  upper  right  corner  of  the 

*  array.   Positive  frequencies  build  from  the  upper  left  corner 

*  diagonally  down  through  the  upper  left  quadrant.    The  other  corners 

*  of  the  array  the  array  are  ALSO  at  the  origin,  so  the  very  center 

*  of  the  array  is  the  HIGHEST  frequency. 
* 

*  In  order  to  use  the  diffraction  filters  we  MUST  take  that  into 

*  account,  that  is  what  the  strange  trigraph  statements  take  care  of 

*  below  */ 

for  (y=0;y<Mcalc;y++) 

{ 

wy  =  (y<M)  ?  y/height  :  (y-Mcalc)/height; 
//  wy  =  (y<M)  ?  y /height  :  (Mcalc-y) /height; 
wy  *  =  wy; 
for(x=0;x<Ncalc;x++){ 

wx=(x<N)  ?  x/width  :   (x-Ncalc)/width; 
//  wx=(x<N)  ?  x/width  :  (Ncalc-x) /width; 

a=lz*(wy+wx*wx)*3. 1415926535; 
/  /  a=  lz*  (wy+  wx*  wx)/F-  PI; 
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ca=cos(a); 
sa=sin(a); 

cc=image[ik]  .re*ca-image[ik]  .im*sa; 
image[ik]  .im=image[ik]  .re*sa+image[ik]  .im*ca; 
image[ik].re=cc; 
ik+  +  ; 
} 
} 
} 

void  OnnralSvnthesizer::MagtoRGB(fTtw_complex  *data) 

{ 

int  i.  j.  ij.  grey,  Ncalc,  Mcalc: 

fftw_real  v.  vmax,  vmin,  vl,  v2; 


Ncalc=2*N; 
Mcalc=2*M; 

hologram = new  RGBImage(N,M); 

vmax=0.0; 
vmin=le6; 
ij=0; 

j  I  Scale  the  result  to  take  maximum  advantage  of  8bit  range 
for(i=0;i<Mcalc*Ncalc;i++){ 

vl=data[i].re; 

v2=data[i].im; 

v=sqrt(vl*vl+v2*v2); 

vmax=MAX(v,  vmax); 
vmin=MIN(v,  vmin); 
} 

vmax=  1  /  (vmax- vmin) ; 


RGBColor  c; 

for(i=0;i<M;i++){ 

for(.j=0;j<N;j  +  +  ){ 
ij=i*Ncalc+j; 

vl=data[ij].im; 
v2=data[ij].re; 

v=(sqrt(vl*vl+v2*v2)-vmin)*vmax; 
data[ij].re=v; 
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grey=(int)floor((fftw_real)(v*255.0)); 
c.B=grey; 
c.R=grey; 
c.G=grey; 

hologram— > plot (j,  i,  c); 

} 
} 


} 


void  OnuralSynthesizer::cleanup() 

{ 

free  (h  plane); 

free(oplane); 

fi'twnd_destroy_plan(fftwFwclPlan); 

fi'twnd_destroy_plan(nrtwBwdPlan); 

} 
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A. 2      Decoder:  Onural  Reconstruction  Algorithm. 

This  is  the  source  code  for  an  Onural  reconstruction  algorithm. 

A. 2.1      OnuralDecoder.hpp 

#ifndef  ONURALDECODER_HPP 
#def  ine  ONURALDECODER_HPP  1 
#include  "Decoder  .hpp'* 
#include  "ProgressUI .hpp" 
#mclude  <fftw.h> 
class  OnuralDecoder:  public  Decoder 

{ 
public: 

OnuralDecoder(); 

virtual  uchar*  decode(uchar  *data,  int  N,  int  M, 

float  W,  float  H,  float  L, 

float  D); 

private: 

int  dJmageN; 
int  d.CalcN; 
int  dJmageM; 
int  d_CalcM; 

float  dJmageW; 
float  dJmageH; 
float  dJmageLambda; 
float  d-decodeDist; 

fft.w .complex  *d_Field; 

uchar  *  d_data; 

void  OnuralDecoder:  :noDClquad(fftw_complex  *field.  fftw_real  average  Value) ; 

fftw_real  OnuralDecoder: :aprinv(fftw_real  zz,  int  Mstage); 

void  OnuralDecoder: :onural(fftw .complex  *hg, 

int  N.  int  M, 
fftw_real  L,  fftw.real  Z, 
fftw_real  X,  fftw_real  Y.int  Mstage); 


fftw.real  Magnitude(fftw .complex  *field); 

void  uchar2field(); 
uchar  *held2uchar(): 
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fftw.real  *clipField(): 


uchar  *OnuralDecoder::mainLoop(); 
ProgressUI  ^progress; 

}; 

#endif 


A. 2. 2      OnuralDecoder.cpp 

#include  "OnuralDecoder .hpp" 
#include  <fftw.h> 
#include  "StatusBox .hpp" 
#include  <FL/F1.H> 
#include  <math.h> 

#def ine  PI  3.14159265359 

#def  ine  TPI  6.28318530718 

#define  MIN(a.  b)  (((a)  <  (b))  ?  (a)  :  (b)) 

#def  ine  MAX(a.  b)  (((a)  >  (b))  ?  (a)  :  (b)) 

OnuralDecoder::OnuralDecoder() 

{ 

d_ImageN=0; 

dJmageM=0; 

d_ImageW=0; 

dJmageH=0; 

d_ImageLambda=0: 

d_decodeDist=0; 

d_data=0; 

progress=new  ProgressUI(); 

} 

uchar  *OnuralDecoder::decode(uchar*  data,  int  N,  int  M,  float  W.  float  H,  float  L. 
float  D) 

{ 

dJmageN=N; 
d_CalcN=2*N; 

dJmageM=M: 
d_CalcM=2*M; 

dJmageW=W; 

dJmageH=H; 

dJmageLambda=L; 

d_decodeDist=D; 

d_data=data; 
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return  mainLoop(); 
} 


/*  Remove  the  DC  component  from  only  the  first  Quadrant,  the  other 
*  three  quadrants  arc  ALREADY  ZERO*/ 

void  OnuralDecoder::noDClquad(fftw .complex  *rield,  fftw_real  averageValue) 

{    . 

int  i,  j,  pos; 

for(i=0;  i<dJmageM;  i++){ 

for(j=0;  j<d_ImageN;  j  +  +  ){ 

pos=i*d_CalcN+j : 

field  [pes]. re  -=  average  Value: 

} 

} 
} 

fftw_real  OnuralDecoder::aprinv(fftw_real  zz,  int  Mstage) 

{ 

fftw_real  cc,ss,powKa; 
int  k,powK: 

cc=cos(zz); 
ss=0; 

if(Mstage>l){ 

ss=2*sin(zz); 
cc+=ss*sin(2*zz); 
ss*  =2*cos(2*zz): 

} 

for  (k=3;k<=Mstage;k++){ 

powK=4«:(k-3);  //  powK=2**(k-l),  fine  for  small  k,  and  FAST 

powKa=(fft,w_real)powK*zz; 

cc+=ss*sin(powKa); 
ss*  =2*cos(powKa): 

} 

return  2*cc/Mstage; 

} 

void  OnuralDecoder::onural(fftw_complex  *hg, 

int  N,  int  M, 

fFtw_real  L,  fftw_real  Z, 

fftw_real  X.  fftw_real  Y.int  Mstage) 
{ 
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int  x,  y.  ml2,  n  12; 
long  ik; 

fftw.real  alpha,  lz.  wx2,  wy2; 

lz=L*Z; 


! 


ml2=M/2; 
nl2=N/2; 

ik=0; 


for  (v=0;y<M;v++)//M 

{ 

wy2  =  (y<ml2)  ?  y/Y  :  (y-M)/Y; 

wy2  *  =  wy2; 

forfx=0;x<N;x++)//7V 

{ 

wx2  =  (x<nl2)  ?  x/X  :  (x-N)/X; 

wx2  *  =  wx2; 

alpha=PI*lz*(wy2+wx2); 

hg[ik].re*  =aprinv(alpha,  Mstage); 
ik++; 

} 

if(  progress— >update(3) )  { 

progress— >-hide(); 
return; 

}: 


} 


fftw_real  OnuralDecoder::Magnitude(fftw_complex  *field) 

{ 

int  i: 

fftw_real  v.sumV,  vl,  v2; 

sumV=0; 

for(i=0;  i  <  d_CalcN*d.CalcM;  i++){ 
vl=field[i].re: 
v2=field[i].im; 
v=sqrt(vl*vl+v2*v2); 
field[i].re=v; 
field[i].im=0; 
sumV+=v: 
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\ 

return  sumV/(d_CalcN  *  cl.CalcM); 

} 

void  OnuralDecoder::uchar2field() 

{ 

int  arraysize=d_CalcN*d_CalcM*sizeof(fftw_complex); 

int  fieldpos,  ucharpos.  x,  y; 
nchar  minval,  maxval; 
fftw_real  scale; 

d_Field=(fftw  .complex  *)calloc(arraysize.  1);/ / calloc  clears  to  zero... 

minval=0; 

maxval=255; 

for(y=0;  y<  dJmageM;  y+  +  ){ 
for(x=0;  x<dJmageN;  x++){ 

ucharpos=3*(x  +  y*d_ImageM); 
minval=MIN(minval,  d_data[ucharpos] ) ; 
maxval=M AX(maxval.  d_data[ucharpos] ) ; 
} 
} 
scale=1.0/(fftw_real)(maxval-minval); 

for(  y=0;  y<  dJmageM;  y++  ){ 
for(  x=0:  x<dJmageN;  x++){ 

fieldpos=x  +  y*d-CalcN; 
\icharpos=3*(x  4-  y*dJmageN); 

d_Field  [fieldpos]. re=((fftw_real)(d_data[ucharpos]-minval))*scale; 
} 
} 

} 

uchar  *OnuralDeooder::field2uchar() 

{ 

int  arraysize=3*dJmageN*d_ImageM*sizeof( uchar); 
int  fieldpos,  ucharpos,  x,  y; 
fftw_real  minval.  maxval,  v,  vre,  vim; 
fftw_real  scale,  val; 

uchar  *  ret  val; 

retval=( uchar*  )calloc(arraysize,  1);// calloc  clears  to  zero... 
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minval=le20; 
maxval=0; 


for(y=0;  y<  dJmageM;  y++  ){ 
for(x=0;  x<d_ImageN;  x++){ 
fieldpos=x  +  y*d_CalcM; 
vre=d_Field[field])os].re; 
vim=d_Field[fieldpos].im; 
v=sqrt(vre*vre+vim*vim); 

minval=MIN(minval,  v); 
maxval=MAX(maxval,  v); 
d_Field  [fieldpos] .  re= v ; 

} 

} 

scale=255.0/(maxval-minval); 

for(y=0;  y<  dJmageM;  y++  ){ 
for(x=0;  x<dJmageN;  x++){ 
fieldpos=x  +  y*d_CalcN; 
ucharpos=3*(x  +  y*d_ImageN); 

val=(uchar)((d_Field  [fieldpos]. re-minval)*scale); 
retval[ucharpos]  =  (uchar)val; 
retval[ucharpos+l]  =  (uchar)val; 
retval[ucharpos+2]  =  (uchar)val; 

} 
} 
return  retval; 


uchar*  OnuralDecoder::mainLoop() 
{ 

fftw-real  DC: 

int  progresscount=0; 

progress— >set.Max(5); 
progress— »update(0) ; 
progress^reset  ( ) ; 
progress— »show(); 

if(progress— »setStatus("Transf  orming  Input.  .  .")){ 
progress->hide(); 
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return  0; 

}; 


ftt.wnd.plan  ForwardFFTPlan,  BackwardFFTPlan; 

uchar2field(); 

progress— mpdate(++progresscount); 
if(progress->setStat,us(  "Making  FFTW  Plans... ")){ 
progress-^hide(); 
return  0; 

}■ 

ForwardFFTPlan=fftw2d_create_plan(d_CalcN,  d_CalcM, 

FFTWJFORWARD, 
FFTW.ESTIM  ATE|FFTW  JN_PLACE) ; 

BackwardFFTPlan=fFt.w2d_create_plan(d_CalcN,  d_CalcM. 

FFTW.BACKWARD, 
FFTW_ESTIMATE|FFTWJN_PLACE); 

DC=4*Magnitude(d_Field);  //The  4  handles  the  doubling  of  array  size 


noDClquad(d_Field,  DC); 

progress— Ripdate(++progresscount): 

if(progress— »setStatus(  "Calculating  Forward  transform..."))! 

progress— »hide(); 

return  0: 

}; 


fi'twnd_one(ForwardFFTPlan,  d_Field,  0); 

d_Field[0].re=0; 

d_Field[0].im=0; 

progress— >update(++progresscount); 
if( progress— >setStatus( "Starting  Onural  reconstruction.  .  .")){ 
progress— >hide(); 
return  0: 

}; 

onural(d_Field.  d_CalcN,  d.CalcM,  dJmageLambda,  d.decodeDist, 
dJmageW,  dJmageH,  5); 
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progress— »update(++progresscount): 
if(progress— »setStatus(  "Calculating  the  Reverse  Transform..."))! 
progress-^hide(); 
return  0; 

}; 


fft.wnd_one(BackwardFFTPlan.  d_Field.  0); 

progress— »update(++progresscount); 
if(  progress —»setStatus(  "Complete"))! 
progress— >  hide  ( ) ; 
return  0; 

}; 

uchar*  decoded=field2uchar(); 

free(d_Field); 

progress— »hide( ) ; 
return  decoded: 
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